Modelling a Zombie Outbreak


Dr Jon Shiach

Senior Lecturer
Department of Computing and Mathematics
Manchester Metropolitan University
Email: j.shiach@mmu.ac.uk
Dr Stephen Lynch

Reader
Department of Computing and Mathematics
Manchester Metropolitan University
Email: s.lynch@mmu.ac.uk
Dr Killian O'Brien

Senior Lecturer
Department of Computing and Mathematics
Manchester Metropolitan University
Email: k.m.obrien@mmu.ac.uk

Introduction

Mathematical modelling is a description of a situation or scenario that is described using mathematical concepts. These models can be very useful as they can allow us to identify possible solutions without first implementing them in the real world. For example, NASA used mathematical models to predict the forces experienced by astronauts to see whether they could survive a trip into space. Here we are going to look at a type of mathematical model called an Agent Based Model (ABM) which can be used to model the behaviour of a group of animals or people. Agent based models consist of a number of autonomous agents which are given a set of rules that govern the agents behaviour in the system. For example, an agent may represent a predator which is given a rule that it will chase the nearest prey. This Jupyter notebook will introduce you to agent based models and how can we write Python programs to model complex behaviour. The programs are written for you but you can make changes to the parameters or even add to the programs if you wish.

Contents

  1. Agent based modelling
  2. Modelling flocking fehaviour
  3. Modelling a zombie outbreak
  4. Modelling the spread of Covid

Jupyter notebooks

Jupyter notebooks are documents that combine text and Python code which allow readable documents such as this one to be created that contain executable code used to perform calculations. To run code in a notebook simply click on the code and then on the run button, or alternatively, press the ctrl + enter keys at the same time. Since we will be using commands to perform calculations and produce our models we need to import some commands to help us to do this. Run the code below to import the commands.

Now we can perform some calculations. The Python code below calculates the length of the hypotenuse of a right angled triangle using the lengths of the two smallest sides and prints the result. Note how the equation $c = \sqrt{a^2 + b^2}$ is entered in a similar way to how we write it on a piece of paper. Can you calculate the length of the hypotenuse when $a = 5$ and $b = 12$


Agent based modelling

An Agent Based Model (ABM) uses individual autonomous agents to which we can give rules that will govern their movement and behaviour. In the models we will use here the agents exist within a finite rectangular region of a given width and height. Each agent is described by its position within the region and velocity. The position is given by a set of co-ordinates $(x,y)$ where $x$ and $y$ are the horizontal and vertical distances from the bottom left-hand corner of the region. The velocity, which governs the direction the agent moves in, is given by the velocity vector (a vector is an arrow with a direction and length) which is represented by $(v_x,v_y)$ where $v_x$ and $v_y$ are the horizontal and vertical speeds of the agent. We can calculate the speed of an agent from the velocity vector using Pythagoras' theorem

$$\textsf{speed} = \sqrt{ v_x ^2 + v_y ^ 2}.$$

This is known as the norm of the velocity vector and is equivalent to the length of the vector. The Python code below defines a function called norm() which calculates the norm of the vector (vx, vy).

To program an agent based model in Python we start by defining an agent object which contains all of the commands that are applied to an agent. The code below defines an agent and gives it a random position and velocity as well as some other attributes which are explained by the comments following the # symbol.

Movement

To model the movement of agents we divide the time which we are simulating into lots of individual time steps. If we know the position and velocity of an agent at one time step we can calculate the position at the next time step using

$$\textsf{new position} = \textsf{current position} + \textsf{velocity vector} \times \textsf{time step}.$$

This is shown in the diagram below where $\mathbf{p}_1$ is the current position of the agent and $\mathbf{p}_2$ and $\mathbf{p}_3$ are the positions in the next two time steps.

In the diagram above we see that the velocity vector does not change so this agent will continue to move along that same path forever. The velocity vector will change if there is some external force applied to the agent, for example, an agent may be attracted towards a particular target. The change in the velocity is known as acceleration which will be determined by any rules that we impose on our agents.

Steering

The steering rules that affect the movement of an agent will each generate a steering vector that are added together to give an acceleration vector

$$\textsf{acceleration vector} = \textsf{steering vector 1} + \textsf{steering vector 2} + \textsf{steering vector 3}$$

An example of a rule imposed on an agent may be that it will travel towards a particular target as shown in the diagram below. The steering vector pointing in the direction of the target is calculated using

$$\textsf{steering vector} = \textsf{target} - \textsf{position}.$$

Since we may be combining several steering rules we may wish to prioritise some over the others. We can do this by limiting the steering vector by a steering factor which is a number that we choose the value of in order to replicate the behaviour we want.

$$\textsf{acceleration vector} = \textsf{steering vector} \times \textsf{steering factor}.$$

The diagram below shows how an agent will move towards a target over several time steps. Note that since the acceleration vector is limited by the steering factor the agents does immediately turn in the direction of the target, rather it gradually turns over a number of time steps.

Keeping the agents within the region

We must ensure that an agent does not pass through the edges of the region and escape. This is shown in the diagram below where the agent is approaching the bottom edge and if no external for is applied it would pass through the edge and out of the region.

To ensure agents remain in the region we apply a steering force if an agent gets too close to an edge which forces it to move away from the edge. A buffer zone is defined that contains the region within some distance of the edge, if an agent is in the buffer zone then a steering force is applied calculated using

\begin{align*} \textsf{steering vector} &= \frac{\textsf{direction vector}}{\textsf{distance from edge}}, \end{align*}

where the direction vector is a vector that points towards the interior of the region. By dividing the direction vector by the distance from the edge the steering force is week when the agent is further away from the edge but becomes progressively stronger as the distance from the decreases. This means that the agent will gradually turn away from the edge.

The code below defines the function avoid_edges() which applies this rule for all four edges of the region defined the by width and height of the region that contains our agents and the size of the buffer zone which is defined as 2 here.


Modelling flocking behaviour

The first model we are going to look at is the modelling of flocking behaviour that can be seen in shoals of fish. Each individual fish in the group reacts to the other fish that surround it and adjusts its movement accordingly. The model Boids which was developed by Craig Reynolds uses three simple rules that are applied to the agents in the model to simulate flocking behaviour, these rules are:

The acceleration forces for all three of these rules is calculated and added to the acceleration vector for an agent.

Neighbouring agents

We want our agents to react to other agents which are the vicinity, for example, if we have an agent that represents a predator it will want to chase after any prey agent and a prey agent will want to flee away from a predator. In order to do this we need to determine which of the other agents can be detected by our agent. The simplest way to do this is to define a detection zone that is a circle surround our agent and all other agents that are within this circle are considered to have been detected by our agent. This could represent a predator using their senses of vision and smell to hunt for their prey. Agents that are within the detection zone are known as neighbouring agents because that are the agents that are nearest to our agent.

The code below defines the function Neighbours() which loops through all of the other agents in the model and calculates the squared distance to a current agent. If this distance is less than the square of the detection radius then it is added to a list of neighbours. The reason why the square of the distance is used is to avoid having to calculate a square root which is slower for a computer to calculate.

Alignment

The alignment rule states that an agent will align itself to the average velocity vector of its neighbouring agents. The average velocity is calculated using

$$\textsf{average velocity} = \frac{\displaystyle\sum \textsf{velocity of neighbour}}{\textsf{number of neighbours}}.$$

The symbol $\sum$ is the Greek character sigma which is equivalent to the Latin character s and is used to denote the sum of a group of numbers. The equation above says that we add up all of the velocity vectors of the neighbouring agents and divide by the number of neighbours to get the average velocity. The steering vector for aligning to the neighbouring agents is then calculated by subtracting the agent velocity from the average velocity of the neighbourng agents

\begin{align*} \textsf{alignment vector} &= \textsf{average velocity} - \textsf{agent velocity}. \end{align*}

The alignment vector is then limited to ensure that its norm is not greater than the alignment factor which controls the extent to which the agent aligns to its neighbours. The code below defines the function alignment() applies the alignment rule to the current agent.

The function limit() used in the alignment() function is defined below. It uses the norm() function to limit a vector v to ensure that its norm is within min_norm and max_norm.

Cohesion

The cohesion rule states that an agent will steer towards the average position of it's neighbouring agents. The steering vector is calculated by subtracting the agent position from the average position of the neigbouring agents

\begin{align*} \textsf{cohesion vector} &= \frac{\displaystyle\sum \textsf{position of neighbour}}{\textsf{number of neighbours}} - \textsf{agent position}. \end{align*}

Like with the alignment vector, the cohesion vector is then limited so that it's norm does not exceed the cohesion factor. The code below defines the function Cohesion() applied the cohesion rule to the current agent.

Separation

The separation rule states that an agent will steer away from agents that are within some small distance of it. For those neighbouring agents that are too close we calculate the steering vector using

\begin{align*} \textsf{separation vector} &= \sum \dfrac{\textsf{position of agent} - \textsf{position of neighbour}}{\textsf{distance to neighbour}}. \end{align*}

By dividing the vector that points away from the neighbouring agent by the distance to that neighbour we scale the separation vector so that the closer a neighbour is, the stronger the seperation force from that neighbour. Once again the separation vector is limited so that it does not exceed the separation factor.

The code below defines a function called separation() applies the separation rule for the current agent.

Applying the agent rules

The code below defines the function agent_rules() which applies the rules that govern the behaviour of the agent. The individual steering vectors for avoiding the edges of the region, alignment, cohesion and separation are calculated and added to the current velocity vector for the agent. The velocity vector is then limited to ensure that it's norm is between the minimum and maximum speed of the agent.

Stepping through time

Now that we have defined the rules for our model we now need to apply these rules for each time step. For each time step we first calculate the rules for each agent in the model and then calculate the new positions of each agent using the acceleration vector that has been calculated when we applied the rules.

The function abm() defined initialises the positions and velocities of the agents in the model and applies the model through time.

Running the model

The following code defines the various parameters that define the model before running the model and generating the animation. Try experimenting with the model by changing these parameters and then clicking on the 'Run' button.


Modelling a zombie outbreak

We have seen with the Boids model that we can use agent based modelling to model a population of animals that exhibit flocking behaviour. We can adapt this model to model the behaviour of human populations. Here we are going to model the (unlikely) event of a zombie outbreak. Zombies have become a popular monster in movies, TV series and computer games such as Night of the Living Dead, The Walking Dead and The Last of Us.

Interpretations of the zombie monster differs so here we will consider a zombie as an reanimated corpse that attempts to feed on humans. If a human is caught and bitten by a zombie they will die and the zombie will feed on the corpse. After a certain amount of time the zombie will cease feeding on the corpse and start looking for another victim and the corpse will eventually reanimate and become a zombie.

To adapt the Boids model we will assign a status to our agents. Initially this is human for most of the agents with a small number of agents being given the status zombie. As the model steps through time the status of a human agent could change to corpse if it has been caught by a zombie and the status of a zombie agent could change to eating as it feeds on a corpse.

The code below defines defines functions that change the status of the agents. Zombies that are eating and corpses are prevented from moving by setting self.move = False and they also given a timer to record the time until they change to the zombie status.

Random path

In the absence of any rules that affect the steering forwards for the agents the agent will just move in a straight line. We don't necessary want this for the zombie model where zombies will wander around so we can make them move along a random path. At each time step, a target is placed at a random position in front of the agent which will cause it to steer towards the target. This is done by calculating the angle that the agent is moving towards relative to a horizontal line and then adding a random angle between $-\pi$ and $\pi$ ($\pi$ is equivalent to 180$^\circ$). Using sine and cosine functions we can then calculate the target position.

The function called RandomPath() defined below alters the agents acceleration vector so that it will move along a random path.

Humans evading zombies

Of course a human agent will want to avoid being eaten by a zombie so we will give the humans a rule that if they detect a zombie, or multiple zombies, they will try to run away in the opposite direction. The code below defines a function called evade() which checks if the agent status is human and if so calculates the acceleration force for evading any neighbouring zombie agents. This is similar to the separation() function used in the Boids model.

Zombie chasing humans

The zombie agents will look for any nearby human agents and chase the nearest one. The code for doing this is shown in the function chase() below. It searches for the nearest human agent and calculates the steering vector that steers the zombie towards the human. If the distance between the zombie and human agents is less than 1 the human has considered to be caught by the zombie so the zombie agent is changed to an eating zombie and the human agent is changed to a corpse.

Applying the agent rules

The agent_rules() function defined below combines the avoid_edges(), random_path(), evade() and chase() rules for our Zombie model.

Running the model

The code below defines the model parameters and runs the model. Here the human agents have been given a maximum speed of 5 whereas the zombie agents have a maximum speed of 3 since common representations of zombies show them as slow lumbering creatures. The zombie creatures in the film 28 Days Later and the computer game The Last of Us move as fast as humans which can easily be modelled by changing the speed parameters. Try experimenting with different parameters to see what affects this has on the model.

Plotting the populations

The animations produced using the abm() function show us which agents are humans and which are zombies at a given time but it isn't easy to see how the number (or population) of each agent types change over time. The function abm_populations() defined below calculates the model and counts the number of each agent type at each time step and stores it in an array called population. After the model has been run it then plots the populations over time.

Note that this function does not generate an animation since that takes too long to compute.


Modelling the spread of Covid

Perhaps a more realistic scenario in which to apply our agent based model is the modelling the spread of the Covid virus. This model is similar to the Zombie model in some respects but without the chase and evade behaviours (we will assume that those people infected with Covid will not want to purposely infect other people).

The Covid virus is like many other viruses such as the common cold where the virus lives within a host. A person infected with the virus may suffer from symptoms such as breathing difficulties and fever which can in rare occasions unfortunately result in the death of the host. After a time the hosts immune system will fight off the virus and the host will recover and become immune from being reinfected by the same virus. During the time when the host is infected with the virus they can infect other people who are close by.

To help prevent the spread of the Covid virus we can implement the following strategies which we can include in our model:

Social distancing can be modelled using the separation() function from the Boids model and quarantining can be modelled by moving an agent out of the region and away from the other agents.

We will start with assigning the majority of our agents with the status susceptible which means that these represent people who have not had Covid so have a possibility of contracting it along with a small number of infected agents. The agents are given a timer attribute which will be used to track has much time has elasped since it was infected. The other types of agents are recovered which are agents that were infected but enough time has elapsed that they have recovered, quarantined which are infected agents that have been placed in quarantine and deceased which are infected agents that have unfortunately died as a result of being infected. This model assumes that recovered agents are immune from reinfection.

Infection

The function infection() defined below used to determine whether an infected agent has passed on the infection to another agent. This is done by checking whether the neighbour is susceptible and whether some random number between 0 and 1 is less than the probability of that the infection is passed. The probability that the infection is passed on is divided by the square of the distance between the infected and susceptible agents so the closer the susceptible agent is the higher the likelihood that they will become infected (and vice-versa). If the infection has been passed to a susceptible agent then we change their status to infected.

Recovery

The function recovery() defined below checks whether an infected agent has recovered. We first subtract the time step from the timer so that it counts down to zero. If the timer is less than zero then the recovery time has elapsed and we change the agent type to that of recovered agent.

Quarantine

The function quarantine() defined below checks whether an infected agent is going into quarantine. We check whether an infected agent needs to enter quarantine by see whether the timer is less than the recovery time $-$ quarantine time.

Fataility

The function fatality() defined below checks whether and infected agent will die within the time step. This is done by checking whether a random number between 0 and 1 is less than the probability that the agent will die within a single time step. If the agent has died its status is changed accordingly.

Agent rules

The function agent_rules() defined below applies the infection(), recovery(), quarantine() fatality() rules for the Covid model as well as the rules avoid_edges() and random_path() for moving the agents. If social distancing is turned on then the separation() rule from the Boids model is used.

Running the model

The code below defines the model parameters and runs the model. Here the social distancing and quarantining have been turned off so the infection should spread quickly. Try experimenting with different parameters to see what affects this has on the model.

Social distancing

To see the affect that social distancing has on the spread of the infection the code below runs the simulation with and without social distancing and produces the population plots.

To see the affects that quarantining has we can produce plots for the model with and without quarantining. What conclusions can you draw from these plots?

Extending the model

The Covid model we have here is quite simplistic and makes the following assumptions about the population.

We could extend our model to take in these additional consideration into account. This is what makes mathematical modelling so useful.


If you would like to further explore Python and fractals you may find the following links useful:

© Dr Jon Shiach 2022